import ipywidgets
%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from tqdm.auto import tqdm
import igraph as ig
import xnet as xn
import importlib
import math
%config InlineBackend.figure_format = 'retina'
from collections import Counter
importlib.reload(xn)

<module 'xnet' from '/N/slate/filsilva/iustudents/xnet.py'>

Student-Student network

Class-Class network

dataSuffix = "2019" #2019 = "", 2020 = "2020"
maxClassSize = 69;
applyFilterMasks = False;
sizeProperty = "CLASS_SIZE_ENROLLED"
# sizeProperty = "CLASS_CAPACITY"

#Loading networks
gstudents = xn.xnet2igraph("Networks/student2student_classes_intersection%s.xnet"%dataSuffix);
# xn.xnet2igraph("Networks/student2student_classes_intersection_withHousing.xnet");
gclasses = xn.xnet2igraph("Networks/classes2classes_students_intersection%s.xnet"%dataSuffix);


#Loading *2* data
students2Classes = [];
classes2Students = [];
# with open("Data/students2Classes.txt","r") as fd:
#   for line in fd:
#     line = line.strip();
#     if(line):
#       entries = line.split("\t");
#       students2Classes.append([int(entry) for entry in entries]);
#     else:
#       students2Classes.append([]);
    
with open("Data/classes2Students%s.txt"%dataSuffix,"r") as fd:
  for line in fd:
    line = line.strip();
    if(line):
      entries = line.split("\t");
      classes2Students.append([int(entry) for entry in entries]);
    else:
      classes2Students.append([]);
gclasses.vs["Students"] = classes2Students;

#Filter classes by size and extra metadata:
maskSize = (np.array(gclasses.vs[sizeProperty])>maxClassSize)

mask = maskSize

if(applyFilterMasks):
  maskMulti = [entry=="YLEC" for entry in gclasses.vs["MULTI_COMPONENTS_FLAG"]]
  maskCLab = [entry.startswith("Y") for entry in gclasses.vs["COMPUTER_LAB_IND"]]
  mask= mask+maskMulti+maskCLab;
  
vertices_toDelete = np.where(mask)[0]
gclasses.delete_vertices(vertices_toDelete);
gclasses.vcount()

#Calculate some properties
gclasses.vs["Strength"] = gclasses.strength(weights="weight")
gclasses.vs["Degree"] = gclasses.degree()
display("calculating closeness...");
gclasses.es["distance"] = 1.0/np.array(gclasses.es["weight"])
# gclasses.vs["Closeness"] = gclasses.closeness(weights="distance")
display("calculating betweenness...");
gclasses.vs["Betweenness"] = gclasses.betweenness(weights="distance",directed=False)
#Random for baseline
gclasses.vs["Random"] = np.random.random(gclasses.vcount())
gclasses.vs["kcore"]= gclasses.coreness()
classes2Students = gclasses.vs["Students"]

'calculating closeness...'

'calculating betweenness...'

40493

%load_ext Cython

The Cython extension is already loaded. To reload it, use:
  %reload_ext Cython

%%cython
#--annotate
cimport numpy as np
import cython
cimport cython
ctypedef np.float64_t DTYPE_t
from libc.stdlib cimport rand, RAND_MAX
from libc.math cimport exp
cdef extern from "stdlib.h":
    double drand48()
    void srand48(long int seedval)

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cdef double randomUniform():
  cdef double r = rand()
  return r / RAND_MAX

def setRandomSeed(long seed):
  srand48(seed);


@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
def updateClass(long[::1] studentStates, double[::1] lastUpdate,long[::1] classesGroupStudents, double gamma,double beta,double beginTime, double endTime):
#   localStates = studentStates[classesGroupStudents];
#   nS = localStates==-1
  cdef long nI = 0;
  cdef long studentIndex;
  cdef double durationSinceUpdate;
  cdef double flip;
  cdef double rRate;
  cdef double duration;
  cdef long oldState;
  for inGroupIndex in range(classesGroupStudents.shape[0]):
    studentIndex = classesGroupStudents[inGroupIndex];
    if(studentStates[studentIndex] == 0):
      nI+=1;
      durationSinceUpdate = endTime - lastUpdate[studentIndex]
      if (durationSinceUpdate>0):
        flip = drand48();
        rRate = 1.0 - exp(-gamma*durationSinceUpdate)
        lastUpdate[studentIndex] = beginTime;
        if flip < rRate:
            studentStates[studentIndex] = 1;
            nI -= 1
  
  for inGroupIndex in range(classesGroupStudents.shape[0]):
    studentIndex = classesGroupStudents[inGroupIndex];
    duration = endTime - beginTime
    oldState = studentStates[studentIndex]
    if (duration>0):
        if (oldState == -1):
            flip = drand48()
            rRate = 1.0 - exp(-beta*duration*nI)
            if flip < rRate:
                lastUpdate[studentIndex] = endTime
                studentStates[studentIndex] = 0
#                 nI += 1 #enable for sync
#             if (oldState == -1):
#                 flip = randomUniform()
#                 rRate = 1.0 - exp(-gamma*duration)
#                 lastUpdate = endTime
#                 if flip < rRate:
#                     studentStates[studentIndex] = 1
#                     nI -= 1

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
def updateRecovery(long[::1] studentStates, double[::1] lastUpdate, double gamma, double currentTime):
#   localStates = studentStates[classesGroupStudents];
#   nS = localStates==-1
  cdef long studentIndex;
  cdef double durationSinceUpdate;
  cdef double flip;
  cdef double rRate;
  for studentIndex in range(studentStates.shape[0]):
    if(studentStates[studentIndex] == 0):
      durationSinceUpdate = currentTime - lastUpdate[studentIndex]
      if (durationSinceUpdate>0):
        flip = drand48();
        rRate = 1.0 - exp(-gamma*durationSinceUpdate)
        lastUpdate[studentIndex] = currentTime;
        if flip < rRate:
            studentStates[studentIndex] = 1;

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
def countSIR(long[::1] studentStates):
#   localStates = studentStates[classesGroupStudents];
#   nS = localStates==-1
  cdef long nS = 0;
  cdef long nI = 0;
  cdef long nR = 0;
  for studentIndex in range(studentStates.shape[0]):
    if(studentStates[studentIndex] == -1):
      nS+=1;
    elif(studentStates[studentIndex] == 0):
      nI+=1;
    elif(studentStates[studentIndex] == 1):
      nR+=1;
  return (nS,nI,nR);

def generateClassesSchedule(p=0.1,propertyName = "Betweenness Centrality"):
  startTimes = gclasses.vs["CLS_START_TIME"]
  endTimes = gclasses.vs["CLS_END_TIME"];
  meetDays = gclasses.vs["CLS_MEETING_DAYS"];
#   classSizes = gclasses.vs["Size"];

  rankingProperty = np.array(gclasses.vs[propertyName]);


  classesByDay = {letter:[] for letter in "MTWRFSN"};
  validClasses = np.where(np.array([d!="nan" for d in meetDays]))[0]
  rankingIndices = rankingProperty[validClasses].argsort()
  maxRank = round((1-p)*len(validClasses))
  selectedIndices = validClasses[rankingIndices][0:maxRank]
  mask = np.ones(gclasses.vcount(), np.bool)
  mask[selectedIndices] = 0
  notSelectedIndices = np.arange(gclasses.vcount())[mask]
  
  gclassesCopy = gclasses.copy();
  gclassesCopy.delete_vertices(notSelectedIndices);
  ggiantClusters = gclassesCopy.clusters()

  numberComponents = len(ggiantClusters)
  sizeLargestComponent = ggiantClusters.giant().vcount()

  for classIndex in selectedIndices:
    if(meetDays[classIndex]!="nan"):
      for weekday in meetDays[classIndex]:
        if(weekday=="D"):
          for wlist in classesByDay:
            classesByDay[wlist].append(classIndex);
        else:
          classesByDay[weekday].append(classIndex);

  #reordering
  studentsInClassesByDay = {};
  startTimesOfClassesByDay = {};
  endTimesOfClassesByDay = {};
  for weekday in classesByDay:
    classesByDay[weekday] = sorted(classesByDay[weekday], key=lambda index: startTimes[index])
    studentsInClassesByDay[weekday] = [np.array(classes2Students[classesIndex],dtype=np.int) for classesIndex in classesByDay[weekday]];
  #   startTimesOfClassesByDay[weekday] = ["%s: %f"%(startTimes[i][-8:-6]+startTimes[i][-5:-3],float(startTimes[i][-8:-6])+float(startTimes[i][-5:-3])/60) for i in classesByDay[weekday]]
    startTimesOfClassesByDay[weekday] = [float(startTimes[i][-8:-6])+float(startTimes[i][-5:-3])/60 for i in classesByDay[weekday]]
    endTimesOfClassesByDay[weekday] = [float(endTimes[i][-8:-6])+float(endTimes[i][-5:-3])/60 for i in classesByDay[weekday]]
  return classesByDay,studentsInClassesByDay,startTimesOfClassesByDay,endTimesOfClassesByDay,numberComponents,sizeLargestComponent

10

weekorder = "MTWRFSN"

def iterateModel(seed=100,gamma = 1.0/6.0/24.0/2.0/(5/7), beta = 0.33/24.0, startProbability = 0.00025, weeksCount = 12, trackBeforeClasses = False):
  setRandomSeed(seed);
  currentDay = 0;
  studentStates = -np.ones(studentsCount,dtype=np.int)
  lastUpdate = np.zeros(studentsCount)
  studentStates[np.random.random(len(studentStates))<startProbability] = 0
  studentHaveClasses = np.zeros(studentsCount,dtype=np.bool)
  for week in weekorder:
    for classesGroupIndex,classesGroupStudents in enumerate(studentsInClassesByDay[week]):
      studentHaveClasses[classesGroupStudents]=True;

  studentStates[~studentHaveClasses]=-2 #Invalid These students do not participate in any classes
  times = [];
  valuesOverTime = [];

  if(not trackBeforeClasses): #First values
    times.append(0.0);
    valuesOverTime.append(countSIR(studentStates));

  for weekIndex in range(weeksCount):
    for week in weekorder:
#       lastClassesStartTime=-10;
      for classesGroupIndex,classesGroupStudents in enumerate(studentsInClassesByDay[week]):
        startTime = currentDay*24+startTimesOfClassesByDay[week][classesGroupIndex];
        endTime = currentDay*24+endTimesOfClassesByDay[week][classesGroupIndex];
        if(trackBeforeClasses):
          times.append(startTime/24.0);
          valuesOverTime.append(countSIR(studentStates));

        updateClass(studentStates, lastUpdate,classesGroupStudents , gamma,beta,startTime, endTime);
      currentDay+=1;
      updateRecovery(studentStates, lastUpdate, gamma, currentDay*24);
      times.append(currentDay);
      valuesOverTime.append(countSIR(studentStates));

  nS,nI,nR = [np.array(entry) for entry in zip(*valuesOverTime)]
  return times,nS,nI,nR;

# from collections import Counter
# studentProperty = gstudents.vs["STUDENT_LEVEL"]

# studentKeys = dict(Counter(studentProperty).most_common(8)).keys()
# studentDistribs = {key:[] for key in studentKeys}
# studentDistribs["Other"] = []

# validStudentProperty = [studentProperty[i] for i,valid in enumerate(studentHaveClasses) if valid];
# currentDistrib = dict(Counter(validStudentProperty))

# remainingStudents = len(validStudentProperty)
# for key in studentKeys:
#   if(key in currentDistrib):
#     keyCount = currentDistrib[key]
#     studentDistribs[key].append(keyCount);
#     remainingStudents-=keyCount;

# studentDistribs["Other"].append(remainingStudents);
# studentDistribs

import multiprocessing
from multiprocessing import Pool
from functools import partial
import random
import sys
from collections import Counter
studentProperty = gstudents.vs["STUDENT_LEVEL"]

gamma = 1.0/6.0/24.0/3.0#/24.0/2.0/(5/7)
beta = 0.33/24.0
startProbability = 0.00025

removedClassesRatio = 0.3
propertyName = "Betweenness"

realizationsCount = 500;


studentsCount = gstudents.vcount();

studentKeys = dict(Counter(studentProperty).most_common(8)).keys()
studentDistribs = {key:[] for key in studentKeys}
studentDistribs["Other"] = []


classesByDay,studentsInClassesByDay,startTimesOfClassesByDay,endTimesOfClassesByDay,numberComponents,sizeLargestComponent = generateClassesSchedule(removedClassesRatio,propertyName)

studentHaveClasses = np.zeros(studentsCount,dtype=np.bool)
for week in weekorder:
  for classesGroupIndex,classesGroupStudents in enumerate(studentsInClassesByDay[week]):
    studentHaveClasses[classesGroupStudents]=True;

validStudentProperty = [studentProperty[i] for i,valid in enumerate(studentHaveClasses) if valid];
currentDistrib = dict(Counter(validStudentProperty))

remainingStudents = len(validStudentProperty)
for key in studentKeys:
  if(key in currentDistrib):
    keyCount = currentDistrib[key]
    studentDistribs[key].append(keyCount);
    remainingStudents-=keyCount;

studentDistribs["Other"].append(remainingStudents);


num_processors = multiprocessing.cpu_count()
p=Pool(processes = num_processors)
allResults = p.map(partial(iterateModel,beta=beta,gamma=gamma,startProbability = startProbability),[np.random.randint(-sys.maxsize-1,sys.maxsize+1) for i in range(realizationsCount)])
p.terminate()

# allResults = []
# for i in tqdm(range(3)):
#   allResults.append(iterateModel(0,seed=np.random.random()));

times,avgS,avgI,avgR = np.average(allResults,axis=0);
_,stdS,stdI,stdR = np.std(allResults,axis=0);
numStudents = np.average(avgS+avgI+avgR);
np.array(allResults)[:,2,:]
peakMagnitudes = np.max(np.array(allResults)[:,2,:],axis=1)
peakPositions = np.argmax(np.array(allResults)[:,2,:],axis=1)

avgPeakMagnitude = np.average(peakMagnitudes)
stdPeakMagnitude = np.std(peakMagnitudes)

totalInfected = (np.array(allResults)[:,2,-1]+np.array(allResults)[:,3,-1])
avgTotalInfected = np.average(totalInfected)
stdTotalInfected = np.std(totalInfected)

avgPeakPosition = np.average(peakPositions)
stdPeakPosition = np.std(peakPositions)

print("Num Students: %d"%numStudents)
print("Num Students2: %d"%np.sum(studentHaveClasses))
print("Total Infected: %f ± %f"%(avgTotalInfected,stdTotalInfected))
print("Peak Magnitude: %f ± %f"%(avgPeakMagnitude,stdPeakMagnitude))
print("Peak Position: %f ± %f"%(avgPeakPosition,stdPeakPosition))
print("Number of comp.: %d"%numberComponents)
print("Size giant comp.: %d"%sizeLargestComponent)
print("Distrib: "+str(studentDistribs))

plt.close("all")
plt.figure()

topS = avgS+stdS
lowS =avgS-stdS
lowS[lowS<0] = 0

topI = avgI+stdI
lowI =avgI-stdI
lowI[lowI<0] = 0

topR = avgR+stdR
lowR =avgR-stdR
lowR[lowR<0] = 0

plt.plot(times,avgS,label="S",color="#848484");
plt.fill_between(times,lowS, topS,alpha=0.25,color="#848484");
plt.plot(times,avgI,label="I",color="#AA403D");
plt.fill_between(times,lowI, topI,alpha=0.25,color="#AA403D");
plt.plot(times,avgR,label="R",color="#7EA0CA");
plt.fill_between(times,lowR, topR,alpha=0.25,color="#7EA0CA");
plt.legend()
plt.xlabel("Days")
plt.show()

Num Students: 36086
Num Students2: 36086
Total Infected: 32573.766000 ± 103.007588
Peak Magnitude: 17371.366000 ± 268.864382
Peak Position: 17.738000 ± 1.540570
Number of comp.: 38
Size giant comp.: 4635
Distrib: {'Senior': [8394], 'Freshman': [7877], 'Sophomore': [7623], 'Junior': [6820], 'Masters': [2530], 'Doctorate': [1581], 'First Year': [295], 'Undergrad Special': [233], 'Other': [733]}

#Over class removal ration


import multiprocessing
from multiprocessing import Pool
from functools import partial
import random
import sys
from collections import Counter

classesRatios = np.linspace(0,1,20,endpoint=False)

gamma = 1.0/6.0/24.0 #/2.0/(5/7)
beta = 0.33/24.0
startProbability = 0.0005

realizationsCount = 100;

studentProperty = gstudents.vs["STUDENT_LEVEL"]

# possibleProperties = ["Strength","Degree","Closeness","Betweenness","Size","kcore","Random"]
possibleProperties = ["Strength","Degree","Betweenness","Size","Random"]

numStudentsArrays={}
avgPeakMagnitudeArrays={}
stdPeakMagnitudeArrays={}
avgTotalInfectedArrays={}
stdTotalInfectedArrays={}
avgPeakPositionArrays={}
stdPeakPositionArrays={}
numberComponentsArrays={}
sizeLargestComponentArrays={}
studentDistribsByProperty={};

pbar = tqdm(total=len(possibleProperties)*len(classesRatios));
for propertyName in tqdm(possibleProperties):
  studentsCount = gstudents.vcount();
  studentKeys = dict(Counter(studentProperty).most_common(8)).keys()
  studentDistribs = {key:[] for key in studentKeys}
  studentDistribs["Other"] = []
  
  numStudentsArray = [];
  avgPeakMagnitudeArray = [];
  stdPeakMagnitudeArray = [];
  avgTotalInfectedArray = [];
  stdTotalInfectedArray = [];
  avgPeakPositionArray = [];
  stdPeakPositionArray = [];
  numberComponentsArray = [];
  sizeLargestComponentArray = [];


  for removedClassesRatio in tqdm(classesRatios):
    pbar.update(1)
    classesByDay,studentsInClassesByDay,startTimesOfClassesByDay,endTimesOfClassesByDay,numberComponents,sizeLargestComponent = generateClassesSchedule(float(removedClassesRatio),propertyName)
  
    studentHaveClasses = np.zeros(studentsCount,dtype=np.bool)
    for week in weekorder:
      for classesGroupIndex,classesGroupStudents in enumerate(studentsInClassesByDay[week]):
        studentHaveClasses[classesGroupStudents]=True;

    validStudentProperty = [studentProperty[i] for i,valid in enumerate(studentHaveClasses) if valid];
    currentDistrib = dict(Counter(validStudentProperty))

    remainingStudents = len(validStudentProperty)
    for key in studentKeys:
      if(key in currentDistrib):
        keyCount = currentDistrib[key]
        studentDistribs[key].append(keyCount);
        remainingStudents-=keyCount;
      else:
        studentDistribs[key].append(0);

    studentDistribs["Other"].append(remainingStudents);


    num_processors = multiprocessing.cpu_count()
    p=Pool(processes = num_processors)
    allResults = p.map(partial(iterateModel,beta=beta,gamma=gamma,startProbability = startProbability),[np.random.randint(-sys.maxsize-1,sys.maxsize+1) for i in range(realizationsCount)])
    p.terminate()

    # allResults = []
    # for i in tqdm(range(3)):
    #   allResults.append(iterateModel(0,seed=np.random.random()));

    times,avgS,avgI,avgR = np.average(allResults,axis=0);
    _,stdS,stdI,stdR = np.std(allResults,axis=0);
    numStudents = np.average(avgS+avgI+avgR);
    
    peakMagnitudes = np.max(np.array(allResults)[:,2,:],axis=1)
    peakPositions = np.argmax(np.array(allResults)[:,2,:],axis=1)
    
    avgPeakMagnitude = np.average(peakMagnitudes)
    stdPeakMagnitude = np.std(peakMagnitudes)
    
    totalInfected = (np.array(allResults)[:,2,-1]+np.array(allResults)[:,3,-1])
    avgTotalInfected = np.average(totalInfected)
    stdTotalInfected = np.std(totalInfected)
    
    avgPeakPosition = np.average(peakPositions)
    stdPeakPosition = np.std(peakPositions)

    numStudentsArray.append(numStudents)
    avgPeakMagnitudeArray.append(avgPeakMagnitude)
    stdPeakMagnitudeArray.append(stdPeakMagnitude)
    avgTotalInfectedArray.append(avgTotalInfected)
    stdTotalInfectedArray.append(stdTotalInfected)
    avgPeakPositionArray.append(avgPeakPosition)
    stdPeakPositionArray.append(stdPeakPosition)
    numberComponentsArray.append(numberComponents)
    sizeLargestComponentArray.append(sizeLargestComponent)

  numStudentsArrays[propertyName] = np.array(numStudentsArray)
  avgPeakMagnitudeArrays[propertyName] = np.array(avgPeakMagnitudeArray)
  stdPeakMagnitudeArrays[propertyName] = np.array(stdPeakMagnitudeArray)
  avgTotalInfectedArrays[propertyName] = np.array(avgTotalInfectedArray)
  stdTotalInfectedArrays[propertyName] = np.array(stdTotalInfectedArray)
  avgPeakPositionArrays[propertyName] = np.array(avgPeakPositionArray)
  stdPeakPositionArrays[propertyName] = np.array(stdPeakPositionArray)
  numberComponentsArrays[propertyName] = np.array(numberComponentsArray)
  sizeLargestComponentArrays[propertyName] = np.array(sizeLargestComponentArray)
  studentDistribsByProperty[propertyName] = studentDistribs

# rows = 4+round(len(possibleProperties)/2)
rows = 2
fig, axs = plt.subplots(rows, 3,figsize=(13,3.0*rows))
title = "Simulation for Beta=%.3f, Gamma=%.3f, Pi=%.2f%%, Max class size = %d, in %s"%(beta*24,gamma*24,startProbability*100,maxClassSize,dataSuffix)
if(applyFilterMasks):
  title+=" (no CL and no LEC in multi. comp.)"
fig.suptitle(title)


plotIndex = 0;
ax = axs.flat[plotIndex];
ax.set_title("Total number of infected")
for propertyName in possibleProperties:
  ax.plot(classesRatios*100,avgTotalInfectedArrays[propertyName],label=propertyName)
  high = avgTotalInfectedArrays[propertyName]+stdTotalInfectedArrays[propertyName]
  low = avgTotalInfectedArrays[propertyName]-stdTotalInfectedArrays[propertyName]
  low[low<0] = 0;
  ax.fill_between(classesRatios*100,low,high,alpha=0.15);
ax.set_xlabel("Removed classes (%)");
ax.set_ylabel("Avg. num. of infected");
ax.set_xlim(0,100)
ax.set_ylim(-500,42000)
plotIndex+=1;

ax = axs.flat[plotIndex];
ax.set_title("Infected peak magnitude")
for propertyName in possibleProperties:
  ax.plot(classesRatios*100,avgPeakMagnitudeArrays[propertyName],label=propertyName)
  high = avgPeakMagnitudeArrays[propertyName]+stdPeakMagnitudeArrays[propertyName]
  low = avgPeakMagnitudeArrays[propertyName]-stdPeakMagnitudeArrays[propertyName]
  low[low<0] = 0;
  ax.fill_between(classesRatios*100,low,high,alpha=0.15);
ax.set_xlabel("Removed classes (%)");
ax.set_ylabel("Avg. num. of infected");
ax.set_xlim(0,100)
ax.set_ylim(-250,20000)
plotIndex+=1;

ax = axs.flat[plotIndex];
ax.set_title("Infected peak position")
for propertyName in possibleProperties:
  ax.plot(classesRatios*100,avgPeakPositionArrays[propertyName],label=propertyName)
  high = avgPeakPositionArrays[propertyName]+stdPeakPositionArrays[propertyName]
  low = avgPeakPositionArrays[propertyName]-stdPeakPositionArrays[propertyName]
  low[low<0] = 0;
  ax.fill_between(classesRatios*100,low,high,alpha=0.15);
ax.set_xlabel("Removed classes (%)");
ax.set_ylabel("Avg. peak position in days");
ax.set_xlim(0,100)
ax.set_ylim(-2,85)
plotIndex+=1;


ax = axs.flat[plotIndex];
ax.set_title("Total number of infected vs Students")
for propertyName in possibleProperties:
  ax.plot(numStudentsArrays[propertyName],avgTotalInfectedArrays[propertyName],label=propertyName)
  high = avgTotalInfectedArrays[propertyName]+stdTotalInfectedArrays[propertyName]
  low = avgTotalInfectedArrays[propertyName]-stdTotalInfectedArrays[propertyName]
  low[low<0] = 0;
  ax.fill_between(numStudentsArrays[propertyName],low,high,alpha=0.15);
ax.set_xlabel("Students with at least 1 class");
ax.set_ylabel("Avg. num. of infected");
ax.set_ylim(-500,40000)
ax.set_xlim(0,42000)
ax.legend()
plotIndex+=1;

ax = axs.flat[plotIndex];
ax.set_title("Infected peak magnitude vs Students")
for propertyName in possibleProperties:
  ax.plot(numStudentsArrays[propertyName],avgPeakMagnitudeArrays[propertyName],label=propertyName)
  high = avgPeakMagnitudeArrays[propertyName]+stdPeakMagnitudeArrays[propertyName]
  low = avgPeakMagnitudeArrays[propertyName]-stdPeakMagnitudeArrays[propertyName]
  low[low<0] = 0;
  ax.fill_between(numStudentsArrays[propertyName],low, high,alpha=0.15);
ax.set_xlabel("Students with at least 1 class");
ax.set_ylabel("Avg. num. of infected");
ax.set_ylim(-250,20000)
ax.set_xlim(0,42000)
plotIndex+=1;

ax = axs.flat[plotIndex];
ax.set_title("Infected peak position vs Students")
for propertyName in possibleProperties:
  ax.plot(numStudentsArrays[propertyName],avgPeakPositionArrays[propertyName],label=propertyName)
  high = avgPeakPositionArrays[propertyName]+stdPeakPositionArrays[propertyName]
  low = avgPeakPositionArrays[propertyName]-stdPeakPositionArrays[propertyName]
  low[low<0] = 0;
  ax.fill_between(numStudentsArrays[propertyName],avgPeakPositionArrays[propertyName]-stdPeakPositionArrays[propertyName], avgPeakPositionArrays[propertyName]+stdPeakPositionArrays[propertyName],alpha=0.15);
ax.set_xlabel("Students with at least 1 class");
ax.set_ylabel("Avg. peak position in days");
ax.set_ylim(-2,85)
ax.set_xlim(0,42000)
plotIndex+=1;

# ax = axs.flat[plotIndex];
# ax.set_title("Num. of connected components")
# for propertyName in possibleProperties:
#   ax.plot(classesRatios*100,numberComponentsArrays[propertyName],label=propertyName)
# ax.set_xlabel("Removed classes (%)");
# ax.set_ylabel("Num. of conn. comp.");
# ax.set_xlim(0,100)
# plotIndex+=1;

# ax = axs.flat[plotIndex];
# ax.set_title("Size of largest component")
# for propertyName in possibleProperties:
#   ax.plot(classesRatios*100,sizeLargestComponentArrays[propertyName],label=propertyName)
# ax.set_xlabel("Removed classes (%)");
# ax.set_ylabel("Classes");
# ax.set_xlim(0,100)
# plotIndex+=1;

# ax = axs.flat[plotIndex];
# for propertyName in possibleProperties:
#   ax.plot(classesRatios*100,numStudentsArrays[propertyName],label=propertyName)
# ax.set_title("Students with at least 1 class")
# ax.set_xlabel("Removed classes (%)");
# ax.set_ylabel("Num. of students");
# ax.set_xlim(0,100)
# plotIndex+=1;

# for propertyName in possibleProperties:
#   ax = axs.flat[plotIndex];
#   ax.set_title("Students by level (%s)"%propertyName)
#   ax.stackplot(classesRatios*100,list(studentDistribsByProperty[propertyName].values()), labels=list(studentDistribsByProperty[propertyName].keys())) 
#   ax.set_xlabel("Removed classes (%)");
#   ax.set_ylabel("Num. of students");
#   ax.set_xlim(0,100)
#   if(plotIndex==5):
#     ax.legend()
#   plotIndex+=1;
  
plt.tight_layout()
fig.subplots_adjust(top=0.87)
plt.show()

nodeStrength = gclasses.strength();
betweennessWeighted = gclasses.vs["Betweenness Centrality"];

bins = np.logspace(np.log10(1), np.log10(np.max(betweennessWeighted)),23)
plt.figure()
plt.plot(nodeStrength,betweennessWeighted,marker="o",ms=2,linestyle = 'None',alpha=0.050)
plt.xscale("log");
plt.yscale("log");
plt.xlabel("Node Strength");
plt.ylabel("Betweenness Centrality");
plt.show()